#include "CParser.h"
#include <fstream>
#include <iostream>
#include <cstdlib>
#include <algorithm>
#include <iomanip>

using namespace std;

//========= local function =================
//compare the clink value
bool CompLink(const CLink *a, const CLink *b)
{
	if(a->i < b->i)		return true;
	else if(a->i == b->i)
		if(a->j < b->j)	return true;
		else			return false;
	else				return false;
}

bool CompPLink(const PLink *a, const PLink *b)
{
	if(a->i < b->i)		return true;
	else if(a->i == b->i)
		if(a->j < b->j)	return true;
		else			return false;
	else return false;
}

bool ReorderPanel(const CPanel *a, const CPanel *b)
{
	if(a->GetCondNum() < b->GetCondNum())		return true;
	else if(a->GetCondNum() == b->GetCondNum())
		if(a->GetCondNum() < b->GetCondNum())	return true;
		else									return false;
	else										return false;
}
//==========================================

//constructor and deconstructor
CParser::CParser()
{
	errM = new double [2];
	errS = new double [2];
}

CParser::~CParser()
{
	delete errM;
	delete errS;
}

void CParser::Start(char **cmd, int size)
{
	Command(cmd, size);		
	RefineScheme();
	OutputImage();
	BuildP();
	SolveCho();
}

//for thesis analysis
void CParser::Start(double lg, double *es, double *em)
{
	Lguard = lg;
	errS = es;
	errM = em;

	RefineScheme();
	BuildP();
	SolveCho();
}

void CParser::Input(char *infile)
{
	char ctemp;
	double coord[12];
	
	panelNum = basisNum = 0;
	
	CPanel *temp_ptr = NULL;		//buffer to vector

	fstream file_ptr;
	file_ptr.open(infile, fstream::in);

	if(!file_ptr.is_open())
	{
		cout<<"Could not open the file!\n"<<endl;
		Help();
		exit(0);
	}

	//input parser
	cout<<"Input file...\n";
	while(file_ptr>>ctemp)
	{
		if(ctemp == 'Q')
		{
			file_ptr>>condNum;

			for(int i=0; i<12; i++)
				file_ptr>>coord[i];

			temp_ptr = new CPanel;
			temp_ptr->InputPanel(0, condNum, coord, panelNum);
			panelptr.push_back(temp_ptr);
			panelNum++;
		}
		else
		{
			//triangle case
		}
	}
	file_ptr.close();

	condNum++;
	InputpanelNum = basisNum = panelNum;
	C = new double [condNum*condNum];
	condVector = new vector<unsigned int> [condNum];

	for(unsigned int n=0; n<InputpanelNum; n++)
		condVector[panelptr[n]->GetCondNum()].push_back(n);
}

//refinement
void CParser::Refine(const int i, const int j, const double *err)
{
//	double Ri = panelptr[i]->GetLongestSideValue();
//	double Rj = panelptr[j]->GetLongestSideValue();

//	if(Ri <= Lguard && Rj <= Lguard)
//		AddCLink(i, j);
//	else
//	{
		if(panelptr[i]->Decision(*panelptr[j], err) && panelptr[j]->Decision(*panelptr[i], err))
			AddCLink(i, j);
		else if(!panelptr[i]->Decision(*panelptr[j], err) && panelptr[j]->Decision(*panelptr[i], err))
		{
			if(!panelptr[i]->isCut())
				panelptr[i]->PanelPartition(panelptr, panelNum, basisNum);
	
			Refine(panelptr[i]->GetLeft(), j, err);
			Refine(panelptr[i]->GetRight(), j, err);
		}
		else if(!panelptr[j]->Decision(*panelptr[i], err) && panelptr[i]->Decision(*panelptr[j], err))
		{
			if(!panelptr[j]->isCut())
				panelptr[j]->PanelPartition(panelptr, panelNum, basisNum);
	
			Refine(panelptr[j]->GetLeft(), i, err);
			Refine(panelptr[j]->GetRight(), i, err);
		}
		else
		{
			if(!panelptr[i]->isCut())
				panelptr[i]->PanelPartition(panelptr, panelNum, basisNum);
		
			if(!panelptr[j]->isCut())
				panelptr[j]->PanelPartition(panelptr, panelNum, basisNum);
		
			Refine(panelptr[i]->GetRight(), panelptr[j]->GetLeft(), err);
			Refine(panelptr[i]->GetLeft(), panelptr[j]->GetLeft(), err);
			Refine(panelptr[i]->GetRight(), panelptr[j]->GetRight(), err);
			Refine(panelptr[i]->GetLeft(), panelptr[j]->GetRight(), err);
		}
//	}
}

//refine scheme
void CParser::RefineScheme()
{
	cout<<"Refining..."<<endl;
	//refinement
	for(unsigned int i=0; i<InputpanelNum; i++)
		for(unsigned int j=i+1; j<InputpanelNum; j++)
		{
			if(panelptr[i]->GetCondNum() != panelptr[j]->GetCondNum())
				Refine(i, j, errM);
		//	else
		//		Refine(i, j, errS);
		}
	
	LeafVector = new vector<unsigned int> [condNum];
	//build basis vector
	for(unsigned int i=0; i<panelNum; i++)
		if(!panelptr[i]->isCut())
			LeafVector[panelptr[i]->GetCondNum()].push_back(i);
	//build self-link
	for(unsigned int n=0; n<condNum; n++)
		for(unsigned int i=0; i<LeafVector[n].size(); i++)
			for(unsigned int j=i; j<LeafVector[n].size(); j++)
				AddCLink(LeafVector[n][i], LeafVector[n][j]);

/*
	for(unsigned int n=0; n<panelNum; n++)
	{
		if(!panelptr[n]->isCut())	AddCLink(n, n);
		if(panelptr[n]->GetNodeType() == 1)	AddCLink(n, n+1);
	}
*/
//	sort(clinkptr.begin(), clinkptr.end(), CompLink);
}

//buile P
void CParser::BuildP()
{
	time_t start, end;
	start = time(NULL);
	cout<<"Build P matrix..."<<endl;
	unsigned int ci, cj, X, Y, entryX, entryY;
	unsigned int i, j, n;
	double cPij;

	for(n=0; n<clinkptr.size(); n++)
	{
		ci = clinkptr[n]->i;
		cj = clinkptr[n]->j;
		cPij = clinkptr[n]->Pij;

		if(panelptr[ci]->GetNodeType() != 2)
		{
			if(panelptr[cj]->GetNodeType() != 2)
			{
				insertP(cPij, ci, cj, false);
			}
			else
			{
				Y = panelptr[cj]->GetParent();
				while(panelptr[Y]->GetNodeType() == 2)
				{
					insertP(-1.0*cPij, ci, panelptr[Y]->GetLeft());		
					Y = panelptr[Y]->GetParent();
				}
				insertP(-1.0*cPij, ci, panelptr[Y]->GetLeft());
				insertP(cPij, ci, Y);
			}
		}
		else	//right
		{
			if(panelptr[cj]->GetNodeType() != 2)
			{	//0:root 1:left
				X = panelptr[ci]->GetParent();
				while(panelptr[X]->GetNodeType() == 2)
				{
					insertP(-1.0*cPij, cj, panelptr[X]->GetLeft());
					X = panelptr[X]->GetParent();
				}

				insertP(-1.0*cPij, cj, panelptr[X]->GetLeft());
				insertP(cPij, cj, X);
			}
			else	//if both panels are not basis
			{
				if(ci != cj)
				{	//not self link
					X = panelptr[ci]->GetParent();
					while(panelptr[X]->GetNodeType() == 2)
					{
						entryX = panelptr[X]->GetLeft();
						Y = panelptr[cj]->GetParent();
						while(panelptr[Y]->GetNodeType() == 2)
						{
							entryY = panelptr[Y]->GetLeft();
							insertP(cPij, entryX, entryY);
							Y = panelptr[Y]->GetParent();
						}
						entryY = panelptr[Y]->GetLeft();
						insertP(cPij, entryX, entryY);
						insertP(-1.0*cPij, entryX, Y);
						X = panelptr[X]->GetParent();
					}
					
					entryX = panelptr[X]->GetLeft();
					Y = panelptr[cj]->GetParent();
					while(panelptr[Y]->GetNodeType() == 2)
					{
						entryY = panelptr[Y]->GetLeft();
						insertP(-1.0*cPij, X, entryY);
						insertP(cPij, entryX, entryY);
						Y = panelptr[Y]->GetParent();
					}

					entryY = panelptr[Y]->GetLeft();
					insertP(cPij, entryX, entryY);
					insertP(-1.0*cPij, entryX, Y);
					insertP(-1.0*cPij, X, entryY);
					insertP(cPij, X, Y);
				}
				else
				{	//self
					X = panelptr[ci]->GetParent();
					while(panelptr[X]->GetNodeType() == 2)
					{
						entryX = panelptr[X]->GetLeft();
						insertP(cPij, entryX, entryX, false);
						Y = panelptr[X]->GetParent();
						while(panelptr[Y]->GetNodeType() == 2)
						{
							entryY = panelptr[Y]->GetLeft();
							insertP(cPij, entryX, entryY, false);
							Y = panelptr[Y]->GetParent();
						}

						entryY = panelptr[Y]->GetLeft();
						insertP(cPij, entryX, entryY, false);
						insertP(-1.0*cPij, entryX, Y, false);
						X = panelptr[X]->GetParent();
					}

					entryX = panelptr[X]->GetLeft();
					insertP(-1.0*cPij, X, entryX, false);
					insertP(cPij, entryX, entryX, false);
					insertP(cPij, X, X, false);
				}
			}
		}
	}
	end = time(NULL);
	cout<<"Time of building P: "<<end-start<<endl;
}

//output image information
void CParser::OutputImage() const
{
	C3DPoint temp;
	ofstream fout_img("Image.txt", ios::out);

	for(unsigned int i=0; i<panelNum; i++)
	{
		fout_img<<"Rectangle"<<endl;
		fout_img<<"panel["<<i<<"] and basis = "<<panelptr[i]->GetBasisIndex();
		fout_img<<" condNum = "<<panelptr[i]->GetCondNum()<<" rootNum = "<<panelptr[i]->GetRoot()<<endl;
		fout_img<<"type = "<<panelptr[i]->GetNodeType()<<endl;
		for(int j=0; j<4; j++)
		{	
			temp = panelptr[i]->GetPanelVertex(j);
			fout_img<<"Point "<<temp.x()<<" "<<temp.y()<<" "<<temp.z()<<endl;
		}
	}
	fout_img.close();
}

//output file for capacitance matrix
void CParser::Output(char *filename)
{
	ofstream fout(filename, ios::out);

	for(int m=0; m<condNum; m++)
	{
		for(int n=0; n<condNum; n++)
		{
			fout.precision (5);
			fout<<setw(8)<<C[m+n*condNum];
		}
		fout<<endl;
	}

	fout.close();
}

//output panel
void CParser::OutputFileForFastCap()
{
	sort(panelptr.begin(), panelptr.end(), ReorderPanel);

	C3DPoint temp;
	ofstream fout("iccap_to_fastcap_test.lst", ios::out);

	fout<<"0 test"<<endl;
	for(unsigned int i=0; i<panelNum; i++)
		if(!panelptr[i]->isCut())
		{
			fout<<"Q "<<panelptr[i]->GetCondNum();
			for(int j=0; j<4; j++)
			{
				temp = panelptr[i]->GetPanelVertex(j);
				fout<<" "<<temp.x()<<" "<<temp.y()<<" "<<temp.z();
			}
			fout<<endl;
		}
	
	fout.close();
}

//Add CLink
void CParser::AddCLink(const unsigned int i, const unsigned int j)
{
	CLink *clinktemp_ptr = new CLink(i, j, panelptr[i]->GetPotential(*(panelptr[j])));
	clinkptr.push_back(clinktemp_ptr);
}

//help
void CParser::Help() const
{
	cout<<"./[a.out] [input_file_name]"<<endl;
	cout<<"[-l] [value] : length guard with input value."<<endl;
	cout<<"[-e] [value] : set eps as input value."<<endl;
	cout<<endl;
}

double CParser::GetPotential(const unsigned int i, const unsigned int j)
{ return panelptr[i]->GetPotential(*(panelptr[j]));}

void CParser::SolveCho()
{
	cout<<"Solve matrix by cholesky decomposition..."<<endl;
	time_t start, end;
	start = time(NULL);
	sort(plinkptr.begin(), plinkptr.end(), CompPLink);

	buff.push_back(plinkptr[0]);
	for(unsigned int n=1; n<plinkptr.size(); n++)
	{
		if(plinkptr[n]->i == plinkptr[n-1]->i && plinkptr[n]->j == plinkptr[n-1]->j)
			buff[buff.size()-1]->ij_value += plinkptr[n]->ij_value;
		else
		{
			if(buff[buff.size()-1]->ij_value == 0)
				buff.pop_back();

			buff.push_back(plinkptr[n]);
		}
	}
	cout<<"element number = "<<2*buff.size()-basisNum<<endl;
	cout<<"basis number = "<<basisNum<<endl;
	cout<<"sparsity = "<<static_cast<double>(2*buff.size()-basisNum)/static_cast<double>(basisNum*basisNum)<<endl;

	choSM P = choSMini(basisNum, buff.size());
	choSM L = NULL;
	choIndex eCount = 0;
	//push items to U matrix
	for(int n=buff.size()-1; n>=0; n--)
		eCount = choPush(P, eCount, basisNum-1-buff[n]->i, basisNum-1-buff[n]->j, buff[n]->ij_value);

	end = time(NULL);
	cout<<"push elements to choSM = "<<end-start<<endl;

	//get Lower matrix of cholesky decomposition
	start = time(NULL);
	L = choleskyFactor(P->vSize, P->xadj, P->adj, P->value);
	end = time(NULL);

	cout<<"time of cholesky decomposition = "<<end-start<<endl;

	//new capacitance matrix
	for(int i=0;i<condNum; i++)
		for(int j=0; j<condNum; j++)
			C[i+j*condNum] = 0;

	//setup V vector and q vector
	double V[basisNum];
	double *q;
	for(int i=0; i<basisNum; i++)
		V[i] = 0.0;

	for(int n=0; n<condNum; n++)
	{
		for(int k=0; k<basisNum; k++)
			V[k] = 0;	//initial
			
		for(int k=0; k<condVector[n].size(); k++)
			V[basisNum-1-condVector[n][k]] = 4*3.1415926*8.85;	

		start = time(NULL);
		q = choSolve(L, V);
		end = time(NULL);
		cout<<"time of solving matrix "<<n<<" = "<<end-start<<endl;

		for(int m=0; m<condNum; m++)
			for(int k=0; k<condVector[n].size(); k++)
				C[m+n*condNum] += q[basisNum-1-condVector[m][k]];
	}
	
	for(int m=0; m<condNum; m++)
	{
		for(int n=0; n<condNum; n++)
		{
			cout.precision (5);
			cout<<setw(8)<<C[m+n*condNum];
		}
		cout<<endl;
	}
}

//insert items in P matrix
void CParser::insertP(double value, unsigned int i, unsigned int j, bool flag)
{
	PLink *plinktmp;
	double val;
	if(i == j)
	{
		val = (!flag) ? value : 2*value;

		plinktmp = new PLink(panelptr[i]->GetBasisIndex(), panelptr[j]->GetBasisIndex(), val);
		plinkptr.push_back(plinktmp);
	}
	else
	{
		if(panelptr[i]->GetBasisIndex() < panelptr[j]->GetBasisIndex())
			plinktmp = new PLink(panelptr[i]->GetBasisIndex(), panelptr[j]->GetBasisIndex(), value);
		else
			plinktmp = new PLink(panelptr[j]->GetBasisIndex(), panelptr[i]->GetBasisIndex(), value);

		plinkptr.push_back(plinktmp);
	}
}

//command line
void CParser::Command(char **cmd, int size)
{
//	Lguard = 0.8;
	errM[0] = 0.01;
	errM[1] = 0.1;
	errS[0] = 0.2;
	errS[1] = 0.3;

	for(int i=0; i<size; i++)
		if(cmd[i][0] == '-')
		{
			if(cmd[i][1] == 'l')
				Lguard = atof(cmd[i+1]);
			if(cmd[i][1] == 'h')
				Help();
			if(cmd[i][1] == 'e')
				if(cmd[i][2] == 'm')
				{
					errM[0] = atof(cmd[i+1]);
					errM[1] = atof(cmd[i+2]);
				}
			if(cmd[i][1] == 'e')
				if(cmd[i][2] == 's')
				{
					errS[0] = atof(cmd[i+1]);
					errS[1] = atof(cmd[i+2]);
				}
		}
}
